#!/usr/bin/python
#coding=utf-8

from __future__ import division, print_function
from ROOT import TCanvas, TH2D, TH1D, TFile, TF1, TGraph
from math import atan

ofile = open('risultati/analisi.dat', 'w')
file_name = 'eventi_calibrati.root'
file = TFile(file_name)
tree = file.Get('eventi_calibrati')

##Istogramma UP
#isto_up = TH1D('isto_up', 'UP; [cm]; [count]', 160, -20, 140)
#canvas_up = TCanvas('canvas_up', 'canvas_up')
#tree.Project('isto_up', 'x_up')
#isto_up.SetStats(0)
#isto_up.Draw()

##Istogramma DOWN
#isto_down = TH1D('isto_down', 'DOWN; [cm]; [count]', 160, -20, 140)
#tree.Project('isto_down', 'x_down')
#canvas_down = TCanvas('canvas_down', 'canvas_down')
#isto_down.SetStats(0)
#isto_down.Draw()

#Istogramma UP vs DOWN
isto_up_vs_down = TH2D('isto_up_vs_down', 'UP vs DOWN; posizione UP [cm]; posizione DOWN [cm]', 160, -20, 140, 160, -20, 140)
tree.Project('isto_up_vs_down', 'x_down:x_up', 'trac==1')
canvas_up_vs_down = TCanvas('canvas_up_vs_down', 'canvas_up_vs_down')
isto_up_vs_down.SetStats(0)
isto_up_vs_down.Draw()
canvas_up_vs_down.SaveAs('./risultati/up_down.eps')

#Istogramma DOWN fissata la posizione sull'UP
H = 100                         #distanza scintillatori
L = 21                          #larghezza scintillatori
x_0s = [10, 60, 110]            #posizione UP
for x_0 in x_0s:
    #Distribuzione teorica con dimensioni geometriche dell'apparato
    z = '(x - {0})'.format(x_0)
    h = 'sqrt(({0}^2 + {1}^2))'.format(z,H)
    w = '({0}/{1})'.format(L,h)
#controllare i segni
    formula = '[0]/{1}^2*({0}*{0}/(1 + {0}^2)^2 + 1/(1 + {0}^2)+ TMath::ATan({0}) - 1)'.format(w,h)
    dist_down_upfixed = TF1('dist_down_upfixed', formula, -20, 140)
    #Distribuzione teorica
    formula2 = '[0]*{0}^2/{1}^2'.format(H, h)
    dist_down_upfixed2 = TF1('dist_down_upfixed2', formula2, -20, 140)
    #Istogramma
    cut_up = 'x_up<({0}+1) && x_up>({0}-1)'.format(x_0)
    isto_down_upfixed = TH1D('isto_down_upfixed', 'DOWN per posizione UP {0} cm; [cm]; [count]'.format(x_0), 40, -20, 140)
    tree.Project('isto_down_upfixed', 'x_down', cut_up)
    canvas_down_upfixed = TCanvas('canvas_down_upfixed', 'canvas_down_upfixed')
    r = isto_down_upfixed.Fit('dist_down_upfixed','rsq', ' ',  0, 120)
    line = 'chi2 down a up fisso a {0} cm: {1}\n'.format(x_0, r.Chi2())
    ofile.write(line)
    isto_down_upfixed.Fit('dist_down_upfixed2','rsq+', ' ',  0, 120)
    isto_down_upfixed.SetStats(1)
    isto_down_upfixed.Draw()
    canvas_down_upfixed.SaveAs('./risultati/down_up_fixed_{0}.eps'.format(x_0))

    cut_down = 'x_down<({0}+1) && x_down>({0}-1)'.format(x_0)
    isto_up_downfixed = TH1D('isto_up_downfixed', 'UP per posizione DOWN {0} cm; [cm]; [count]'.format(x_0), 40, -20, 140)
    tree.Project('isto_up_downfixed', 'x_up', cut_down)
    canvas_up_downfixed = TCanvas('canvas_up_downfixed', 'canvas_up_downfixed')
    r = isto_up_downfixed.Fit('dist_down_upfixed','rsq', ' ',  0, 120)
    line = 'chi2 up a down fisso a {0} cm: {1}\n'.format(x_0, r.Chi2())
    ofile.write(line)
    #isto_down_upfixed.Fit('dist_down_upfixed2','rsq+', ' ',  0, 120)
    isto_up_downfixed.SetStats(1)
    isto_up_downfixed.Draw()
    canvas_up_downfixed.SaveAs('./risultati/up_down_fixed_{0}.eps'.format(x_0))

file.Close()


##Distribuzione angolare
angles = [0, 30, 60, 90]
time = 1500                         #tempo di acquisizione [s]
file_gr_int = open('gr_int.dat', 'w')
for angle in angles:
    file = TFile('./dati/{0}_gr.root'.format(angle))
    tree = file.Get('pjmca')
    events = tree.GetEntries()
    intensity = events/time         #[raggi/s]
    line = '{0}\t{1}\n'.format(angle, intensity)
    file_gr_int.write(line)
    file.Close()
file_gr_int.close()
graph_gr_int = TGraph('gr_int.dat')
canvas_gr_int = TCanvas('canvas_gr_int', 'canvas_gr_int')
graph_gr_int.SetMarkerStyle(2)
delta_ang_rad = atan(L/H)
dist_ang_teo = '[0]*({0} + 1/4*(sin(2*({1} + {0}))- sin(2*({1} - {0}))))'.format(delta_ang_rad, 'x*TMath::DegToRad()')
dist_teo = TF1('dist_teo', dist_ang_teo, -10, 100)
r = graph_gr_int.Fit(dist_teo, 'sq', '', -10, 100)
line = 'chi2 dist ang: {0}\n'.format(r.Chi2())
ofile.write(line)
graph_gr_int.SetTitle('Distribuzione angolare;angolo[deg];frequeza dei raggi [Hz]')
graph_gr_int.Draw('ap')
canvas_gr_int.SaveAs('./risultati/gr_int.eps')

ofile.close()
